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Abstract 

The long-time behavior of certain fast-decaying infinite temper- 
ature correlation functions on one-, two- and three-dimensional lat- 
tices of classical spins with various kinds of nearest-neighbor inter- 
actions is studied numerically, and evidence is presented that the 
functional form of this behavior is either simple exponential or expo- 
nential multiplied by cosine. Due to the fast characteristic timescale 
of the long-time decay, such a universality cannot be explained on 
the basis of conventional Markovian assumptions. It is suggested 
that this behavior is related to the chaotic properties of the spin 
dynamics. 
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1 Introduction 



We perform a numerical study of the fong-time behavior of infinite temper- 
ature correlation fmictions defined on an infinite lattice of classical spins 
as: 

Git) = ( S-,{t) E„cos(q . r,.„) 5,^(0) ), (1) 

where Sj^ is the /ith or z) spin component on the fcth lattice site; rfc„ 
is the translation vector between the fcth and the ?ith sites; and q is a 
wave vector commensurate with the lattice periodicity. We consider three 
types of lattices: a simple one-dimensional chain, a two-dimensional square 
lattice, and a three-dimensional cubic lattice. In each case, the dynami- 
cal evolution of the system is driven by the nearest-neighbor interaction 
represented by the Hamiltonian 

H = ^[JjjS'^S'J^ + JySfS^^ + JzS^S^, (2) 

where are coupling constants. With such a Hamiltonian, the timescale of 
the individual spin motion referred to below as the "mean free time" can be 
given by the time t = [\NS'^{Jx^ + Jy^ + J^)Y^^'^-, where N is the number 
of the nearest neighbors (twice the number of the lattice dimensions) . 

In the context of inelastic neutron scattering, correlation functions |^ 
are called "intermediate structure factors" Pj. If g = 0, Eq.QJ can also 
represent the free induction decay in nuclear magnetic resonance (NMR)j2]. 

In this work, we provide extensive numerical evidence that the generic 
long-time behavior of G(t) has one of the following two functional forms: 
either 

G{t) 2. e-«*, (3) 

or 

G(i) ~ e-«*cos(77i-f 0), (4) 

where the constants ^ and 77 are of the order of 1/t. 

It is important to realize that, if the above functional form of the long- 
time behavior is, indeed, generic (i.e. independent of the specific details 
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of interaction) , then this property is very likely related to the randomness 
generated by the spin dynamics. At the same time, the problem cannot be 
reduced to the Markovian paradigm of "a slow variable interacting with a 
fast equilibrating background" — the characteristic timescale r in Eas. H3l4|) 
is not "slow". It is, in fact, the fastest natural timescale of the problem. 
Therefore, whatever is the ultimate explanation of that behavior, it will 
certainly be a step beyond the standard theory of Brownian-type motion. 

Our interest in the long-time behavior of the correlation functions 1^ 
was originally motivated by two isolated pieces of evidence supporting the 
oscillatory behavior in quantum (spin 1/2) systems: (i) experiments 
on NMR free induction decay in CaF2 [3] and (ii) the results of numerical 
diagonalization of spin 1/2 chains 0]. In the both cases, quantities analo- 
gous to the one defined by Eq.JQ) have been measured or computed, and 
the results look very similar to the plots shown in the left column of Fig. 1. 

We came to recognize the importance of a detailed study of the classical 
limit, when, in an attempt to explain the long-time relaxation in quantum 
spin systems, we developed a theory that turned out to be simultaneously 
applicable to classical spins. That theory is presented in a different paper|5] 
which has been written simultaneously with the present one. The present 
paper is mainly numerical: it is not intended to be a brief exposition of 
Ref.|ni. Below, we only provide the summary of the results from Ref.jS]. 

The theory developed in Ref. 151 describes long-time relaxation as a cor- 
related diffusion in finite volumes. In the classical case, those finite volumes 
correspond to the spherical surfaces on which the tips of classical spin vec- 
tors move, while in the quantum case, the finite volumes originate from a 
more sophisticated construction in Hilbert space. The overall structure of 
such a treatment has noticeable parallels with the theory of PoUicott-Ruelle 
resonances in classical chaotic systems [HJ [7]. A definite prediction from the 
correlated diffusion description is that the functional form of the long-time 
relaxation should be given by Eas. (|3l4| l. 

The important part of the above theory is not the diffusion descrip- 
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tion itself but the reason why it is apphcable, given the "non-Markovian" 
relaxation timescale. The theory is based on the fairly strong conjecture 
that, for a broad class of many-body systems, a formal extension of the 
Brownian-like description applies to the the long-time behavior of the en- 
semble average quantities, even when the problem exhibits no separation 
of the timescales between the slow and the fast motions. 

Before proceeding with the description of the simulations, it should be 
mentioned that, for the classical spin systems at infinite temperature, the 
long-time behavior of the g-dependent correlation functions decaying 
on the timescale of r has never been addressed. The closest to this subject 
was the work of de Alcantara Bonfim and ReiterlSj, who considered the 
Heisenberg spin chain and focused on the long-time behavior of correla- 
tion functions ^ with small q. Those correlation hmctions, however, are 
not typical for our purposes, because, as a consequence of the total spin 
conservation, they decay on the timescale, which is much longer than the 
characteristic timescale of one-spin motion. In that situation, the hypoth- 
esis of spin diffusion[5] would lead to the prediction of nearly exponential 
decay with the decay constant proportional to . The results of de Alcan- 
tara Bonfim and Rcitcr did not cover the range of values, which would be 
sufficient to confirm or rule out the exponential character of the long-time 
decay. However, those results (in line with others^n],^^) indicated that, 
if, the spin diffusion regime exists for classical spin chains, the approach to 
that regime is anomalously slow. 

The present work includes one example of the Heisenberg interaction, 
just to show that this case does not appear to be special with respect to 
the long-time property (|3I4|I . 

2 Simulations 

Our computational strategy was similar to that of Miiller 10 . Namely, we 
did not deal with very large systems but, instead, performed an ensemble 
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averaging over a large number of finite, but not too small, lattices having 
periodic boundary conditions. The finite size effects were then controlled 
by varying the size of the lattice. 

For a given lattice size, many computational runs have thus been per- 
formed. Each of them started from completely random initial conditions 
(corresponding to the infinite temperature) and generated the evolution of 
the system over a time interval two orders of magnitude longer than the 
mean free time r (see Table 1 for specific numbers). The correlation func- 
tions were then obtained by averaging the data within each run and over 
different runs. 

The following algorithm has been used in order to simulate the evolution 
of the system. 

At each time step, the spins were advanced sequentially in such a way 
that, if the spin number k interacted with the spin number n, and the A;th 
spin was advanced first, then the new coordinates of the nth spin were 
computed based on the local field created by the already advanced fcth 
spin. 

The procedure for advancing a given (fcth) spin to the next point along 
the discrete time grid consisted of two steps. 

Step 1: The coordinates of the A;th spin were changed by 6Sk according 
to the straightforward discretization of the equations of motion, i.e. 

6Sk = [Sk X hk] 5t, (5) 

where hk was the local field equal to J2^'^'[JxS^ex + JyS^Cy + JzS^Cz] 
(sum over the nearest neighbors of the A;th spin); and St was the discretiza- 
tion time step. 

Step 2: The higher order errors that changed the length of the spin 
vector were eliminated. This was done by contracting the spin component 
perpendicular to the local field, so that it took the absolute value it had 
before Step 1. 

The whole manipulation could not change the spin projection parallel 
to the local field and, therefore, the energy of interaction of that spin with 
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its neighbors. Since the next spin was advanced in the newly updated local 
field, the energy of the whole system was conserved exactly during the entire 
integration. Apart from insuring meaningful behavior of the computed 
trajectories, the exact conservation of energy substantially improved the 
convergence of the algorithm with respect to the limit 6t — > 0. 

In our simulations, the discretization time steps (given in Table 1) were 
admitted as sufficient when their further reduction appeared to have no 
effect on the computed correlation functions. We also checked that the 
averaging over a larger number of much shorter runs led to results consistent 
with the longer runs we used. 

The time lengths of the runs indicated in Table 1 were chosen to opti- 
mize the resulting efficiency of averaging: too short runs (of the order of 
the time length of the computed correlation function) did not make many 
independent contributions to the correlation functions, while too long runs 
did not improve the quality of the averaging proportionally to their length. 

3 Results 

The results of our simulations for different dimensions, interaction con- 
stants, and wave numbers are presented in Fig. 1 jl2j. together with the 
long-time theoretical fits based on either Eq.© or Since we aimed at 
demonstrating the exponential character of the long-time behavior, it was 
natural to use a logarithmic scale for G{t). However, because in the half of 
the cases, G{t) was also oscillating, we chose to show the logarithmic plots 
for the absolute value of G{t), which explains the cusps in the left column 
of plots in Fig. 1. Those cusps correspond to the points where G{t) crosses 
zero. The reason that the cusp minima do not reach — cx) is that a discrete 
grid was used. 

The selection of parameters for the simulations was subject to cer- 
tain practical constraints: The computed correlation functions could be 
considered reliable only within quite a limited range along both the t- 
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and Log(G(i))- axes. Therefore, we avoided the correlation functions that 
reached too-small values too fast, or, on the contrary, decayed too slowly. 

From our experience, the finite size effects were least pronounced for 
correlation functions with q — 0. For this reason, q was chosen to be zero 
in six of eight examples presented in Fig. ^ 

Since the long-time behavior of the correlation functions was only marginally 
accessible with our computational resources, we present the results in sub- 
stantial detail, thus making clear the uncertainties associated with insuf- 
ficient ensemble averaging and finite size effects. Each of the correlation 
functions presented in Fig. ^ was computed four times: two statistically 
independent averaging results for each of two different lattice sizes. "Two 
statistically independent averaging results" means that, in each case, the 
same number of sample runs was performed but two different sets of ran- 
dom numbers were used for setting the initial orientations of spins. Each 
frame in that figure thus contains a superposition of four plots. The time 
interval where these plots do not deviate from each other can be consid- 
ered as representing the limit of infinite lattice size with sufficient ensemble 
averaging. 

The finite size effects are not evident in any of the examples shown in 
Fig. 1 — in every case, the plots representing different simulation outcomes 
do not deviate from each other before the statistical fluctuations for each 
of the two lattice sizes become apparent. 

Our experience indicates that further improvement in the accuracy of 
the computed correlation functions would simultaneously require a finer dis- 
cretization, much more extensive ensemble averaging and, probably, larger 
system sizes, i.e. much greater computational effort. 

Summarizing the evidence, we observe that, with the marginal excep- 
tion of Fig. 1(d), in every other case presented, there is an interval, covering 
at least one decade of the values of G{t), where the simulation results agree 
with the long-time fits © or Q). We would also like to point out that, while 
in all cases the exponential behavior becomes pronounced quite early, the 
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long-time exponents in Figs. l(e,g) describe almost the entire correlation 
functions. 

Thus our numerical results lend strong support to the idea expressed in 
Ref. |S] , that a discrete spectrum of well-separated exponents describes the 
long-time behavior of the correlation functions considered, with the slowest 
of those exponents responsible for the asymptotic functional form given by 
Eq.lPl) or It is also very likely, though slightly less reliable, that in 
the examples presented, our simulations revealed the slowest exponents. 
The reservation here is for the possibility that even slower exponents could 
enter the long-time expansion of G(t) with anomalously small coefficients. 

4 Conclusions 

In conclusion, we have presented a numerical evidence that the long-time 
behavior of the correlation functions considered is exponential with or with- 
out the oscillatory component. Leaving the details to Ref.^, here we just 
mention that our best hope for the theoretical explanation of that behavior 
is associated with the strong chaotic properties of the spin dynamics. Those 
properties are likely to be quite generic, i.e. present in other systems. 
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TABLE CAPTIONS: 
Table 1: 

Simulation parameters and the long-time fits corresponding to the plots 
presented in Fig. ^ The numbers in the right four columns characterize 
each of the four data sets superimposed in the corresponding frame. 
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Frame Lattice dimensions Number Time length Discreti- 
label of runs of each run, zation 

Size 1 Size 2 [{JSy^] time step, 



(a) 


19 


15 


(b) 


19 


15 


(c) 


40 


24 


(d) 


10 X 10 


7x7 


(e) 


10 X 10 


7x7 


(f) 


12 X 8 


8x4 


(g) 


5x5x5 


4x4x3 


(h) 


5x5x5 


4x4x3 


Frame 


Long-time fit 


label 


(time units 





500000 320 0.02 

500000 320 0.02 

160000 320 0.02 

60000 100 0.0025 

80000 170 0.01 

60000 100 0.0025 

80000 170 0.005 

80000 170 0.005 



(a) 2.25 cxp(-0.585 t) cos(1.93 t + 0.31) 

(b) 0.23 exp(-0.357 t) 

(c) 0.30 exp(-0.588 t) cos(1.46 t - 1.19) 

(d) 0.40 exp(-0.645 t) 

(e) 1.20 exp(-1.031 t) cos(3.06 t - 0.82) 

(f) 0.88 exp(-0.488 t) 

(g) 1.00 exp(-1.299 t) cos(3.70 t - 0.82) 

(h) 1.70 exp(-0.426 t) 

Table 1: 
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FIGURE CAPTIONS: 



Figure 1: 

Correlation functions G{t) of the form ^ for ID chain, 2D square 
lattice, and 3D cubic lattice. The interaction coefficients and the wave 
numbers (in the units of inverse lattice spacing) are indicated above the 
plots. The main frame of each figure shows the logarithmic scale of the 
absolute value of G{t)/G{0), while the inset frame shows the direct plots of 
G{t)/G{0) (with neither the logarithm nor the absolute value being taken). 
Within each frame, the simulation results are presented by almost indis- 
tinguishable superposition of data for two lattice sizes, and each size is 
represented by two statistically independent averaging results; therefore, 
two solid lines for the larger size (Size 1) and two dash-dotted lines for the 
smaller size (Size 2). The spread of the four lines indicates the computa- 
tional uncertainty — it becomes visible only in the lower right corner of 
each frame. The dashed lines in each figure are the long-time theoretical 
fits of form or The numbers relevant to each data set are given in 
Table 1. 
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(b) 1D: Jx= 1.2J; Jy= -0.3J; Jz= J. 



(C) 1D: Jx= J; Jy= J; Jz= J. 



(e) 2D: Jx=0;Jy=-J;Jz=J. 



(g) 3D: 4=0; Jy= -J; Jz= J. 







(d) 2D: Jx=1.2J;Jy=-0.3J;Jz=J. 

q=(0,0) 




(f) 2D: Jx= 1.2J; Jy= -0.2J; Jz= J. 
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(h) 3D: Jx=-0.5J;Jy=0.5J;Jz=J. 




t, [units of (JS)"^ ] 
Figure 1: 
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